library(survey)
library(foreign)
library(ggplot2)
library(cowplot)
library(dplyr)
library(egg)
library(ggpubr)
library(gridGraphics)

## see wave2_CLEANING.R for data cleaning
load("~/Dropbox/Public_Conf_Mil/Data_and_code/wave2_data/clean_dataw2.RData")

## Create weighted survey design object
w2_design <-
  svydesign(
    id = ~ 1,
    weights = ~ weight2,
    data = df
  )

#### FIGURE 10.1 - DIFF IN MEANS TREATMENT EFFECTS #### 

iran_df <- data.frame(out = c(NA, svyttest(Q21_d ~ A2_support, w2_design)$estimate, NA,
                              svyttest(Q21_d ~ A2_oppose, w2_design)$estimate, NA),
                      ci_lo = c(NA, svyttest(Q21_d ~ A2_support, w2_design)$conf.int[1], NA,
                                svyttest(Q21_d ~ A2_oppose, w2_design)$conf.int[1], NA),
                      ci_hi = c(NA, svyttest(Q21_d ~ A2_support, w2_design)$conf.int[2], NA,
                                svyttest(Q21_d ~ A2_oppose, w2_design)$conf.int[2], NA),
                      cond = c(1:5))

trns_df <- data.frame(out = c(NA, svyttest(Q25_d ~ A3_support, w2_design)$estimate, NA,
                              svyttest(Q25_d ~ A3_oppose, w2_design)$estimate, NA,
                              svyttest(Q25_d ~ A3_divide, w2_design)$estimate, NA),
                      ci_lo = c(NA, svyttest(Q25_d ~ A3_support, w2_design)$conf.int[1], NA,
                                svyttest(Q25_d ~ A3_oppose, w2_design)$conf.int[1], NA,
                                svyttest(Q25_d ~ A3_divide, w2_design)$conf.int[1], NA),
                      ci_hi = c(NA, svyttest(Q25_d ~ A3_support, w2_design)$conf.int[2], NA,
                                svyttest(Q25_d ~ A3_oppose, w2_design)$conf.int[2], NA,
                                svyttest(Q25_d ~ A3_divide, w2_design)$conf.int[2], NA),
                      cond = c(1:7))

iran_plot <- ggplot(iran_df) + 
  geom_pointrange(aes(x = cond, y = out, ymin = ci_lo, ymax = ci_hi), fatten = 2, size = 0.5) +
  geom_text(aes(x = cond, y = out, label = sprintf("%0.2f", round(out, digits = 2))), 
            size = 3,
            vjust = 2.5) + xlab("Treatment Condition") + 
  scale_x_continuous(breaks = c(2,4),
                     labels = c("Support", "Oppose")) +
  ylab("Difference in Means from Control Condition") + ylim(-0.2, 0.2) +
  geom_hline(yintercept = 0, colour = "red", lty = 1) + theme_bw() +
  ggtitle("Iran Strikes")
iran_plot

trns_plot <- ggplot(trns_df) + 
  geom_pointrange(aes(x = cond, y = out, ymin = ci_lo, ymax = ci_hi), fatten = 2, size = 0.5) +
  geom_text(aes(x = cond, y = out, label = sprintf("%0.2f", round(out, digits = 2))), 
            size = 3,
            vjust = 2.5) + xlab("Treatment Condition") + 
  scale_x_continuous(breaks = c(2,4,6),
                     labels = c("Support", "Oppose", "Divided")) +
  ylab("Difference in Means from Control Condition") + ylim(-0.2, 0.2) +
  geom_hline(yintercept = 0, colour = "red", lty = 1) + theme_bw() +
  ggtitle("Transgender Ban")
trns_plot


fig10_1 <- ggarrange(iran_plot + theme(axis.text.y = element_text(size = 6),
                                       plot.title = element_text(size = 11),
                                       axis.title.y = element_text(size = 7),
                                       axis.title.x = element_text(size = 9)), 
                     trns_plot + theme(axis.text.y = element_text(size = 6),
                                       plot.title = element_text(size = 11),
                                       axis.title.y = element_text(size = 7),
                                       axis.title.x = element_text(size = 9)),
                     nrow = 1, ncol = 2)
fig10_1
ggsave("~/Dropbox/Public_Conf_Mil/Chapter 10 Figures/fig10_1_treatments.jpeg", 
       plot = fig10_1, dpi = 700)

#### FIGURE 10.2 - DIFF IN MEANS TREATMENT EFFECTS, HANDGUNS IN SCHOOLS ####

guns_df <- data.frame(out = c(svyttest(Q27_d ~ A4_more, w2_design)$estimate,
                              svyttest(Q27_d ~ A4_less, w2_design)$estimate,
                              svyttest(Q27_d ~ A4_divd, w2_design)$estimate),
                      ci_lo = c(svyttest(Q27_d ~ A4_more, w2_design)$conf.int[1],
                                svyttest(Q27_d ~ A4_less, w2_design)$conf.int[1],
                                svyttest(Q27_d ~ A4_divd, w2_design)$conf.int[1]),
                      ci_hi = c(svyttest(Q27_d ~ A4_more, w2_design)$conf.int[2],
                                svyttest(Q27_d ~ A4_less, w2_design)$conf.int[2],
                                svyttest(Q27_d ~ A4_divd, w2_design)$conf.int[2]),
                      cond = c(1:3))

fig_guns <- ggplot(guns_df) + 
  geom_pointrange(aes(x = cond, y = out, ymin = ci_lo, ymax = ci_hi), fatten = 2, size = 0.5) +
  geom_text(aes(x = cond, y = out, label = sprintf("%0.2f", round(out, digits = 2))), 
            size = 3,
            vjust = 2.5) + xlab("Treatment Condition") + 
  scale_x_continuous(breaks = c(1:3),
                     labels = c("More Safe", "Less Safe",
                                "Divided")) +
  ylab("Difference in Means from Control Condition") + ylim(-0.1, 0.1) +
  geom_hline(yintercept = 0, colour = "red", lty = 1) + theme_bw()
fig_guns
ggsave("~/Dropbox/Public_Conf_Mil/Chapter 10 Figures/fig10_2_handguns.jpeg", 
       plot = fig_guns, dpi = 700)

#### FIGURE 10.3 - TREATMENT EFFECTS BY CONFIDENCE ####

iran_ef <- svyby(~ Q21_d, ~ P_ASSIGN2 + Q12_d, w2_design, svymean, na.rm = TRUE)
trns_ef <- svyby(~ Q25_d, ~ P_ASSIGN3 + Q12_d, w2_design, svymean, na.rm = TRUE)

iran_tr <- data.frame(out = c(iran_ef[c(1,4,5),3] * 100, NA, iran_ef[c(6,9,10),3] * 100),
                      ci_lo = c(iran_ef[c(1,4,5),3] * 100 - 1.96*iran_ef[c(1,4,5),4] * 100, NA,
                                iran_ef[c(6,9,10),3] * 100 - 1.96*iran_ef[c(6,9,10),4] * 100),
                      ci_hi = c(iran_ef[c(1,4,5),3] * 100 + 1.96*iran_ef[c(1,4,5),4] * 100, NA,
                                iran_ef[c(6,9,10),3] * 100 + 1.96*iran_ef[c(6,9,10),4] * 100),
                      cond = c(1:7))

trns_tr <- data.frame(out = c(trns_ef[c(1,3,2,4),3] * 100, NA, trns_ef[c(5,7,6,8),3] * 100),
                      ci_lo = c(trns_ef[c(1,3,2,4),3] * 100 - 1.96*trns_ef[c(1,3,2,4),4] * 100, NA,
                                trns_ef[c(5,7,6,8),3] * 100 - 1.96*trns_ef[c(5,7,6,8),4] * 100),
                      ci_hi = c(trns_ef[c(1,3,2,4),3] * 100 + 1.96*trns_ef[c(1,3,2,4),4] * 100, NA,
                                trns_ef[c(5,7,6,8),3] * 100 + 1.96*trns_ef[c(5,7,6,8),4] * 100),
                      cond = c(1:9))

fig_iran <- ggplot(iran_tr) + theme_bw() +
  geom_col(aes(x = cond, y = out), width = 0.75, alpha = 0.7) +
  geom_errorbar(aes(x = cond, ymin = ci_lo, ymax = ci_hi), size = 0.5, alpha = 0.6, width = 0.4) +
  scale_x_continuous(breaks = c(1:7),
                     labels = c("  Control", "Support", "Oppose", "",
                                "Control", "Support", "Oppose")) +
  geom_text(aes(x = cond, y = out, label = sprintf("%0.1f", round(out, digits = 1))),
            size = 2.5,
            vjust = -0.2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  ylab("Respondent Agreement (%)") + ylim(0,50) + 
  xlab("Not Confident                            Confident") +
  ggtitle("Iran Strikes") 
fig_iran

fig_trns <- ggplot(trns_tr) + theme_bw() +
  geom_col(aes(x = cond, y = out), width = 0.75, alpha = 0.7) +
  geom_errorbar(aes(x = cond, ymin = ci_lo, ymax = ci_hi), size = 0.5, alpha = 0.6, width = 0.4) +
  scale_x_continuous(breaks = c(1:9),
                     labels = c("  Control", "Support", "Oppose", "Divided", "",
                                "Control", "Support", "Oppose", "Divided")) +
  geom_text(aes(x = cond, y = out, label = sprintf("%0.1f", round(out, digits = 1))),
            size = 2.5,
            vjust = -0.2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  ylab("Respondent Agreement (%)") + ylim(0,50) + 
  xlab("Not Confident                            Confident") +
  ggtitle("Transgender Ban") 
fig_trns

fig10_3_a <- ggarrange(fig_iran + theme(axis.text.y = element_text(size = 6),
                                        plot.title = element_text(size = 11),
                                        axis.title.y = element_text(size = 9),
                                        axis.text.x = element_text(size = 8),
                                        axis.title.x = element_text(size = 9)), 
                       fig_trns + theme(axis.text.y = element_text(size = 6),
                                        plot.title = element_text(size = 11),
                                        axis.title.y = element_text(size = 9),
                                        axis.text.x = element_text(size = 8),
                                        axis.title.x = element_text(size = 9)),
                       nrow = 1, ncol = 2)
fig10_3_a
ggsave("~/Dropbox/Public_Conf_Mil/Chapter 10 Figures/fig10_3_barchart.jpeg", 
       plot = fig10_3_a, dpi = 700)

#### FIGURE 10.3 - DIFF IN MEANS TREATMENT EFFECTS BY CONFIDENCE ####

iran_df2 <- data.frame(out = c(svyttest(Q21_d ~ A2_support_n, w2_design)$estimate,
                               svyttest(Q21_d ~ A2_oppose_n, w2_design)$estimate, NA,
                               svyttest(Q21_d ~ A2_support_c, w2_design)$estimate,
                               svyttest(Q21_d ~ A2_oppose_c, w2_design)$estimate),
                       ci_lo = c(svyttest(Q21_d ~ A2_support_n, w2_design)$conf.int[1],
                                 svyttest(Q21_d ~ A2_oppose_n, w2_design)$conf.int[1], NA,
                                 svyttest(Q21_d ~ A2_support_c, w2_design)$conf.int[1],
                                 svyttest(Q21_d ~ A2_oppose_c, w2_design)$conf.int[1]),
                       ci_hi = c(svyttest(Q21_d ~ A2_support_n, w2_design)$conf.int[2],
                                 svyttest(Q21_d ~ A2_oppose_n, w2_design)$conf.int[2], NA,
                                 svyttest(Q21_d ~ A2_support_c, w2_design)$conf.int[2],
                                 svyttest(Q21_d ~ A2_oppose_c, w2_design)$conf.int[2]),
                       cond = c(1:5))

trns_df2 <- data.frame(out = c(svyttest(Q25_d ~ A3_support_n, w2_design)$estimate,
                               svyttest(Q25_d ~ A3_oppose_n, w2_design)$estimate, NA,
                               svyttest(Q25_d ~ A3_support_c, w2_design)$estimate,
                               svyttest(Q25_d ~ A3_oppose_c, w2_design)$estimate),
                       ci_lo = c(svyttest(Q25_d ~ A3_support_n, w2_design)$conf.int[1],
                                 svyttest(Q25_d ~ A3_oppose_n, w2_design)$conf.int[1], NA,
                                 svyttest(Q25_d ~ A3_support_c, w2_design)$conf.int[1],
                                 svyttest(Q25_d ~ A3_oppose_c, w2_design)$conf.int[1]),
                       ci_hi = c(svyttest(Q25_d ~ A3_support_n, w2_design)$conf.int[2],
                                 svyttest(Q25_d ~ A3_oppose_n, w2_design)$conf.int[2], NA,
                                 svyttest(Q25_d ~ A3_support_c, w2_design)$conf.int[2],
                                 svyttest(Q25_d ~ A3_oppose_c, w2_design)$conf.int[2]),
                       cond = c(1:5))

iran_plot2 <- ggplot(iran_df2) + 
  geom_hline(yintercept = 0, colour = "red", lty = 1) + theme_bw() +
  geom_pointrange(aes(x = cond, y = out, ymin = ci_lo, ymax = ci_hi), fatten = 2, size = 0.5) +
  geom_text(aes(x = cond, y = out, label = sprintf("%0.2f", round(out, digits = 2))), 
            size = 3,
            vjust = 2.5) + xlab("Treatment Condition") + 
  scale_x_continuous(breaks = c(1,2,4,5),
                     labels = c("Support", "Oppose", "Support", "Oppose")) +
  ylab("Difference in Means from Control Condition") + ylim(-0.2, 0.2) +
  xlab("Confident                                            Not Confident") + 
  ggtitle("Iran Strikes")
iran_plot2

trns_plot2 <- ggplot(trns_df2) + 
  geom_hline(yintercept = 0, colour = "red", lty = 1) + theme_bw() +
  geom_pointrange(aes(x = cond, y = out, ymin = ci_lo, ymax = ci_hi), fatten = 2, size = 0.5) +
  geom_text(aes(x = cond, y = out, label = sprintf("%0.2f", round(out, digits = 2))), 
            size = 3,
            vjust = 2.5) + xlab("Treatment Condition") + 
  scale_x_continuous(breaks = c(1,2,4,5),
                     labels = c("Support", "Oppose", "Support", "Oppose")) +
  ylab("Difference in Means from Control Condition") + ylim(-0.2, 0.2) +
  xlab("Confident                                            Not Confident") + 
  ggtitle("Transgender Ban")
trns_plot2

fig10_3 <- ggarrange(iran_plot2 + theme(axis.text.y = element_text(size = 6),
                                        plot.title = element_text(size = 11),
                                        axis.title.y = element_text(size = 7),
                                        axis.text.x = element_text(size = 6),
                                        axis.title.x = element_text(size = 8)), 
                     trns_plot2 + theme(axis.text.y = element_text(size = 6),
                                        plot.title = element_text(size = 11),
                                        axis.title.y = element_text(size = 7),
                                        axis.text.x = element_text(size = 6),
                                        axis.title.x = element_text(size = 8)),
                     nrow = 1, ncol = 2)
fig10_3
ggsave("~/Dropbox/Public_Conf_Mil/Chapter 10 Figures/fig10_3_treatbyconf.jpeg", 
       plot = fig10_3, dpi = 1000)

#### FIGURE 10.4 - TRANSGENDER BY BY CONFIDENCE, PARTY ####

trns_full <- svyby(~ Q25_d, ~ P_ASSIGN3 + party + Q12_d, w2_design, svymean, na.rm = TRUE)

trns_tr2 <- data.frame(dem_n = c(trns_full[c(1,3,2,4),4]*100),
                       dem_n_lo = c(trns_full[c(1,3,2,4),4]*100 - 
                                      1.96*trns_full[c(1,3,2,4),5]*100),
                       dem_n_hi = c(trns_full[c(1,3,2,4),4]*100 + 
                                      1.96*trns_full[c(1,3,2,4),5]*100),
                       dem_c = c(trns_full[c(13,15,14,16),4]*100),
                       dem_c_lo = c(trns_full[c(13,15,14,16),4]*100 - 
                                      1.96*trns_full[c(13,15,14,16),5]*100),
                       dem_c_hi = c(trns_full[c(13,15,14,16),4]*100 + 
                                      1.96*trns_full[c(13,15,14,16),5]*100),
                       rep_n = c(trns_full[c(9,11,10,12),4]*100),
                       rep_n_lo = c(trns_full[c(9,11,10,12),4]*100 -
                                      1.96*trns_full[c(9,11,10,12),5]*100),
                       rep_n_hi = c(trns_full[c(9,11,10,12),4]*100 +
                                      1.96*trns_full[c(9,11,10,12),5]*100),
                       rep_c = c(trns_full[c(21,23,22,24),4]*100),
                       rep_c_lo = c(trns_full[c(21,23,22,24),4]*100 -
                                      1.96*trns_full[c(21,23,22,24),5]*100),
                       rep_c_hi = c(trns_full[c(21,23,22,24),4]*100 +
                                      1.96*trns_full[c(21,23,22,24),5]*100),
                       dem_n_cond = c(1:4),
                       dem_c_cond = c(12:15),
                       rep_n_cond = c(6:9),
                       rep_c_cond = c(17:20))
trns_tr2$dem_n_lo[3] <- 0

fig_trns2 <- ggplot(trns_tr2) + theme_bw() +
  geom_col(aes(x = dem_n_cond, y = dem_n), width = 0.75, alpha = 0.7, color = "blue", fill = "blue") +
  geom_col(aes(x = dem_c_cond, y = dem_c), width = 0.75, alpha = 0.7, color = "blue", fill = "blue") +
  geom_col(aes(x = rep_n_cond, y = rep_n), width = 0.75, alpha = 0.7, color = "red", fill = "red") +
  geom_col(aes(x = rep_c_cond, y = rep_c), width = 0.75, alpha = 0.7, color = "red", fill = "red") +
  geom_errorbar(aes(x = dem_n_cond, ymin = dem_n_lo, ymax = dem_n_hi), size = 0.5, 
                alpha = 0.6, width = 0.4) +
  geom_errorbar(aes(x = dem_c_cond, ymin = dem_c_lo, ymax = dem_c_hi), size = 0.5, 
                alpha = 0.6, width = 0.4) +
  geom_errorbar(aes(x = rep_n_cond, ymin = rep_n_lo, ymax = rep_n_hi), size = 0.5, 
                alpha = 0.6, width = 0.4) +
  geom_errorbar(aes(x = rep_c_cond, ymin = rep_c_lo, ymax = rep_c_hi), size = 0.5, 
                alpha = 0.6, width = 0.4) +
  geom_text(aes(x = dem_n_cond, y = dem_n, label = sprintf("%0.1f", round(dem_n, digits = 1))),
            size = 2.5, vjust = -0.2) +
  geom_text(aes(x = dem_c_cond, y = dem_c, label = sprintf("%0.1f", round(dem_c, digits = 1))),
            size = 2.5, vjust = -0.2) +
  geom_text(aes(x = rep_n_cond, y = rep_n, label = sprintf("%0.1f", round(rep_n, digits = 1))),
            size = 2.5, vjust = -0.2) +
  geom_text(aes(x = rep_c_cond, y = rep_c, label = sprintf("%0.1f", round(rep_c, digits = 1))),
            size = 2.5, vjust = -0.2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  scale_x_continuous(breaks = c(1:20),
                     labels = c("  Control", "Support", "Oppose", "Divided", "",
                                "Control", "Support", "Oppose", "Divided", "", "",
                                "Control", "Support", "Oppose", "Divided", "",
                                "Control", "Support", "Oppose", "Divided")) +
  ylab("Respondent Agreement (%)") + ylim(0,80) + 
  xlab("Not Confident                                                        Confident")
fig_trns2
ggsave("~/Dropbox/Public_Conf_Mil/Chapter 10 Figures/fig10_4_barchart.jpeg", 
       plot = fig_trns2, dpi = 700)

#### FIGURE 10.4 - DIFF IN MEANS TRANSGENDER BY CONFIDENCE, PARTY ####

trns_df3 <- data.frame(rep = c(svyttest(Q25_d ~ A3_support_n_rep, w2_design)$estimate,
                               svyttest(Q25_d ~ A3_oppose_n_rep, w2_design)$estimate,
                               svyttest(Q25_d ~ A3_support_c_rep, w2_design)$estimate,
                               svyttest(Q25_d ~ A3_oppose_c_rep, w2_design)$estimate),
                       rep_lo = c(svyttest(Q25_d ~ A3_support_n_rep, w2_design)$conf.int[1],
                                  svyttest(Q25_d ~ A3_oppose_n_rep, w2_design)$conf.int[1],
                                  svyttest(Q25_d ~ A3_support_c_rep, w2_design)$conf.int[1],
                                  svyttest(Q25_d ~ A3_oppose_c_rep, w2_design)$conf.int[1]),
                       rep_hi = c(svyttest(Q25_d ~ A3_support_n_rep, w2_design)$conf.int[2],
                                  svyttest(Q25_d ~ A3_oppose_n_rep, w2_design)$conf.int[2],
                                  svyttest(Q25_d ~ A3_support_c_rep, w2_design)$conf.int[2],
                                  svyttest(Q25_d ~ A3_oppose_c_rep, w2_design)$conf.int[2]),
                       dem = c(svyttest(Q25_d ~ A3_support_n_dem, w2_design)$estimate,
                               svyttest(Q25_d ~ A3_oppose_n_dem, w2_design)$estimate,
                               svyttest(Q25_d ~ A3_support_c_dem, w2_design)$estimate,
                               svyttest(Q25_d ~ A3_oppose_c_dem, w2_design)$estimate),
                       dem_lo = c(svyttest(Q25_d ~ A3_support_n_dem, w2_design)$conf.int[1],
                                  svyttest(Q25_d ~ A3_oppose_n_dem, w2_design)$conf.int[1],
                                  svyttest(Q25_d ~ A3_support_c_dem, w2_design)$conf.int[1],
                                  svyttest(Q25_d ~ A3_oppose_c_dem, w2_design)$conf.int[1]),
                       dem_hi = c(svyttest(Q25_d ~ A3_support_n_dem, w2_design)$conf.int[2],
                                  svyttest(Q25_d ~ A3_oppose_n_dem, w2_design)$conf.int[2],
                                  svyttest(Q25_d ~ A3_support_c_dem, w2_design)$conf.int[2],
                                  svyttest(Q25_d ~ A3_oppose_c_dem, w2_design)$conf.int[2]),
                       rcond = c(1,4,8,11),
                       dcond = c(2,5,9,12))

trns_plot3 <- ggplot(trns_df3) + 
  geom_pointrange(aes(x = rcond, y = rep, ymin = rep_lo, ymax = rep_hi), color = "red") +
  geom_pointrange(aes(x = dcond, y = dem, ymin = dem_lo, ymax = dem_hi), color = "blue") +
  geom_text(aes(x = rcond, y = rep, label = sprintf("%0.2f", round(rep, digits = 2))), size = 3,
            vjust = 2) +
  geom_text(aes(x = dcond, y = dem, label = sprintf("%0.2f", round(dem, digits = 2))), size = 3,
            vjust = 2) + 
  xlab("Not Confident                                                               Confident") + 
  scale_x_continuous(breaks = c(1.5, 4.5, 8.5, 11.5),
                     labels = c("Support", "Oppose", "Support", "Oppose")) +
  ylab("Difference in Means from Control Condition") +
  geom_hline(yintercept = 0, colour = "black", lty = 1) +
  theme_bw()
trns_plot3
ggsave("~/Dropbox/Public_Conf_Mil/Chapter 10 Figures/fig10_4_ban_conf_party.jpeg", 
       plot = trns_plot3, dpi = 1000)

#### FIGURE 10.5 - PLANNING, ADVICE, INTEGRATION ####

plan_c <- svyby(~ Q24A_d, ~ P_LEAD + party, w2_design, svymean, na.rm = TRUE) # treatments
advc_c <- svyby(~ Q24D_d, ~ P_LEAD + party, w2_design, svymean, na.rm = TRUE)
intg_c <- svyby(~ Q24E_d, ~ P_LEAD + party, w2_design, svymean, na.rm = TRUE)

plan_m <- svyby(~ Q24B_d, ~ party, w2_design, svymean, na.rm = TRUE) # military (all respondents)
advc_m <- svyby(~ Q24C_d, ~ party, w2_design, svymean, na.rm = TRUE)
intg_m <- svyby(~ Q24F_d, ~ party, w2_design, svymean, na.rm = TRUE)

plan_df <- data.frame(out = c(plan_c[1:3,3] * 100, plan_m[1,2] * 100, NA, 
                              plan_c[4:6,3] * 100, plan_m[2,2] * 100, NA,
                              plan_c[7:9,3] * 100, plan_m[3,2] * 100),
                      ci_lo = c(plan_c[1:3,3] * 100 - 1.96*plan_c[1:3,4] * 100, 
                                plan_m[1,2] * 100 - 1.96*plan_m[1,3] * 100, NA, 
                                plan_c[4:6,3] * 100 - 1.96*plan_c[4:6,4] * 100, 
                                plan_m[2,2] * 100 - 1.96*plan_m[2,3] * 100, NA,
                                plan_c[7:9,3] * 100 - 1.96*plan_c[7:9,4] * 100, 
                                plan_m[3,2] * 100 - 1.96*plan_m[3,3] * 100),
                      ci_hi = c(plan_c[1:3,3] * 100 + 1.96*plan_c[1:3,4] * 100, 
                                plan_m[1,2] * 100 + 1.96*plan_m[1,3] * 100, NA, 
                                plan_c[4:6,3] * 100 + 1.96*plan_c[4:6,4] * 100, 
                                plan_m[2,2] * 100 + 1.96*plan_m[2,3] * 100, NA,
                                plan_c[7:9,3] * 100 + 1.96*plan_c[7:9,4] * 100, 
                                plan_m[3,2] * 100 + 1.96*plan_m[3,3] * 100),
                      cond = c(1:14))

advc_df <- data.frame(out = c(advc_c[1:3,3] * 100, advc_m[1,2] * 100, NA, 
                              advc_c[4:6,3] * 100, advc_m[2,2] * 100, NA,
                              advc_c[7:9,3] * 100, advc_m[3,2] * 100),
                      ci_lo = c(advc_c[1:3,3] * 100 - 1.96*advc_c[1:3,4] * 100, 
                                advc_m[1,2] * 100 - 1.96*advc_m[1,3] * 100, NA, 
                                advc_c[4:6,3] * 100 - 1.96*advc_c[4:6,4] * 100, 
                                advc_m[2,2] * 100 - 1.96*advc_m[2,3] * 100, NA,
                                advc_c[7:9,3] * 100 - 1.96*advc_c[7:9,4] * 100, 
                                advc_m[3,2] * 100 - 1.96*advc_m[3,3] * 100),
                      ci_hi = c(advc_c[1:3,3] * 100 + 1.96*advc_c[1:3,4] * 100, 
                                advc_m[1,2] * 100 + 1.96*advc_m[1,3] * 100, NA, 
                                advc_c[4:6,3] * 100 + 1.96*advc_c[4:6,4] * 100, 
                                advc_m[2,2] * 100 + 1.96*advc_m[2,3] * 100, NA,
                                advc_c[7:9,3] * 100 + 1.96*advc_c[7:9,4] * 100, 
                                advc_m[3,2] * 100 + 1.96*advc_m[3,3] * 100),
                      cond = c(1:14))

intg_df <- data.frame(out = c(intg_c[1:3,3] * 100, intg_m[1,2] * 100, NA, 
                              intg_c[4:6,3] * 100, intg_m[2,2] * 100, NA,
                              intg_c[7:9,3] * 100, intg_m[3,2] * 100),
                      ci_lo = c(intg_c[1:3,3] * 100 - 1.96*intg_c[1:3,4] * 100, 
                                intg_m[1,2] * 100 - 1.96*intg_m[1,3] * 100, NA, 
                                intg_c[4:6,3] * 100 - 1.96*intg_c[4:6,4] * 100, 
                                intg_m[2,2] * 100 - 1.96*intg_m[2,3] * 100, NA,
                                intg_c[7:9,3] * 100 - 1.96*intg_c[7:9,4] * 100, 
                                intg_m[3,2] * 100 - 1.96*intg_m[3,3] * 100),
                      ci_hi = c(intg_c[1:3,3] * 100 + 1.96*intg_c[1:3,4] * 100, 
                                intg_m[1,2] * 100 + 1.96*intg_m[1,3] * 100, NA, 
                                intg_c[4:6,3] * 100 + 1.96*intg_c[4:6,4] * 100, 
                                intg_m[2,2] * 100 + 1.96*intg_m[2,3] * 100, NA,
                                intg_c[7:9,3] * 100 + 1.96*intg_c[7:9,4] * 100, 
                                intg_m[3,2] * 100 + 1.96*intg_m[3,3] * 100),
                      cond = c(1:14))

fig_plan <- ggplot(plan_df) + theme_bw() +
  geom_col(aes(x = cond, y = out), width = 0.75, alpha = 0.7) +
  geom_errorbar(aes(x = cond, ymin = ci_lo, ymax = ci_hi), size = 0.5, alpha = 0.6, width = 0.4) +
  scale_x_continuous(breaks = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14),
                     labels = c("Civilian", "Democratic", "    Republican", "Military", "",
                                "Civilian", "Democratic", "    Republican", "Military", "",
                                "Civilian", "Democratic", "    Republican", "Military")) +
  geom_text(aes(x = cond, y = out, label = sprintf("%0.1f", round(out, digits = 1))),
            size = 2.5,
            vjust = -0.2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  ylab("Respondent Agreement (%)") + ylim(0,65) +
  xlab("Democrat Respondents                     Independent Respondents                    Republican Respondents") +
  ggtitle("Political leaders had a good plan/Military leaders implemented plan well") 
fig_plan

fig_advc <- ggplot(advc_df) + theme_bw() +
  geom_col(aes(x = cond, y = out), width = 0.75, alpha = 0.7) +
  geom_errorbar(aes(x = cond, ymin = ci_lo, ymax = ci_hi), size = 0.5, alpha = 0.6, width = 0.4) +
  scale_x_continuous(breaks = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14),
                     labels = c("Civilian", "Democratic", "    Republican", "Military", "",
                                "Civilian", "Democratic", "    Republican", "Military", "",
                                "Civilian", "Democratic", "    Republican", "Military")) +
  geom_text(aes(x = cond, y = out, label = sprintf("%0.1f", round(out, digits = 1))),
            size = 2.5,
            vjust = -0.2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  ylab("Respondent Agreement (%)") + ylim(0,65) + 
  xlab("Democrat Respondents                    Independent Respondents                     Republican Respondents") +
  ggtitle("Political leaders listened to military advice/Military leaders gave good advice") 
fig_advc

fig_intg <- ggplot(intg_df) + theme_bw() +
  geom_col(aes(x = cond, y = out), width = 0.75, alpha = 0.7) +
  geom_errorbar(aes(x = cond, ymin = ci_lo, ymax = ci_hi), size = 0.5, alpha = 0.6, width = 0.4) +
  scale_x_continuous(breaks = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14),
                     labels = c("Civilian", "Democratic", "    Republican", "Military", "",
                                "Civilian", "Democratic", "    Republican", "Military", "",
                                "Civilian", "Democratic", "    Republican", "Military")) +
  geom_text(aes(x = cond, y = out, label = sprintf("%0.1f", round(out, digits = 1))),
            size = 2.5,
            vjust = -0.2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  ylab("Respondent Agreement (%)") + ylim(0,65) + 
  xlab("Democrat Respondents                    Independent Respondents                      Republican Respondents") +
  ggtitle("Political/Military leaders did a good job integrating military, diplomatic, and economic efforts")
fig_intg

fig10_5 <- ggarrange(fig_plan + theme(axis.text.y = element_text(size = 6),
                                        plot.title = element_text(size = 11),
                                        axis.title.y = element_text(size = 7),
                                        axis.text.x = element_text(size = 8),
                                        axis.title.x = element_text(size = 9)), 
                     fig_advc + theme(axis.text.y = element_text(size = 6),
                                        plot.title = element_text(size = 11),
                                        axis.title.y = element_text(size = 7),
                                        axis.text.x = element_text(size = 8),
                                        axis.title.x = element_text(size = 9),
                                      plot.margin = unit(c(3,3,3,3), "lines")),
                     fig_intg + theme(axis.text.y = element_text(size = 6),
                                        plot.title = element_text(size = 11),
                                        axis.title.y = element_text(size = 7),
                                        axis.text.x = element_text(size = 8),
                                        axis.title.x = element_text(size = 9)), 
                     nrow = 3, ncol = 1)
fig10_5
ggsave("~/Dropbox/Public_Conf_Mil/Chapter 10 Figures/fig10_5_planning.jpeg", 
       plot = fig10_5, dpi = 700)
